Transcriptome-based identification and validation of reference genes for corm growth stages, different tissues, and drought stress in Taro (Colocasia esculenta)

Taro is a widely utilized starch resource plant. It is essential to quantify the expression levels of functional genes associated with taro growth using real-time quantitative polymerase chain reaction (RT-qPCR). However, to obtain reliable RT-qPCR results, appropriate reference genes (RGs) are required for data normalization. In this study, we screened seven novel candidate RGs using transcriptome datasets from taro, encompassing data from growth corms and various tissues. The expression stability of these seven new RGs, along with the commonly used RGs Actin, EF1-α, and β-tubulin, was assessed using Delta Ct, BestKeeper, geNorm, and NormFinder algorithms. Furthermore, we conducted a comprehensive analysis using the RefFinder program and validated the results using the target gene, CeAGPL1. The findings revealed that ACY-1 and PIA2 were the optimal multiple RGs for normalization during corm growth, while COX10 and Armc8 were suitable for samples including various types of tissues. Furthermore, we found three RGs, Armc8, COX10 and CCX4L, were the optimal RGs for drought stress. This study assessed the suitability of RGs in taro for the first time. The identified RGs provide valuable resources for studying corm growth, diverse tissues, and drought stress. This study contributes to the advancement of our understanding of the underlying mechanisms that govern the growth of taro. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-024-05199-x.


Introduction
Taro (Colocasia esculenta) is a starchy root crop that has been cultivated and consumed by humans for thousands of years [1].It is a tropical and subtropical plant that belongs to the Araceae family and is widely grown in many countries around the world, including Asia, Africa, the Caribbean, and the Pacific Islands [2,3].Taro is known by various names in different regions, such as "dasheen" in the Caribbean [3], "kalo" in Hawaii [4], "arvi" in India [5], and "eddo" in Japan [6].Taro is a staple food for many indigenous communities and has played a significant role in their cultural and dietary practices.It is a versatile crop that can be used in a variety of culinary applications, such as boiling, frying, baking, and mashing.Taro leaves are edible and are commonly used in traditional dishes.In addition to its culinary importance, taro has been used in traditional medicine because of its potential health benefits [7].Apart from its cultural and dietary significance, taro is also economically important as a cash crop.It is a source of income for many smallscale farmers in developing countries, and the global trade in taro and taro-derived products has been steadily increasing.Taro is gaining popularity in modern cuisines because of its unique flavor, nutritional value, and potential as a gluten-free alternative to wheat-based products [8].
Taro corm, which is the edible part, contains high levels of carbohydrates, fats, crude fiber, vitamin C, thiamin, riboflavin, and niacin [9].With a rich carbohydrate content and energy value of 135 kcal/100 g, taro corms have almost double the carbohydrate content of potatoes, and their protein content is 11% higher than that of yams, cassava, and sweet potatoes [10].The growth of taro corm occurs in three phases.Given the nutrient-rich profile of taro corms, it is important to understand their genetic basis for the corm growth, particularly the regulation of pathways related to carbohydrate accumulation [11,12].
Despite its importance, research on taro has been limited compared to that on other major crops.However, with advancements in molecular biology and genomics, there has been increasing interest in studying taro at the molecular level to understand the genetics, physiology, and molecular pathways associated with important traits, such as yield, disease resistance, and nutritional content.One crucial aspect of molecular research is the validation of reference genes (RGs), which are used as internal controls for normalizing gene expression data in quantitative gene expression studies, such as reverse transcription quantitative polymerase chain reaction (RT-qPCR).Actin was used as a control gene to determine the transcript levels of a metal-responsive gene, pCeMT, in C. esculenta [13].Lekshmi et al. used actin as a RG to quantify the expression of resistant gene analogue in taro with qPCR [14].In additional, Wang et al. normalized the expression of multiple genes during cold storage using the actin7 gene [15].In a transcriptome study of taro corm growth, actin was selected as an RG to verify the transcriptome.Obviously, researchers tended to use actin as RG in taro studies.Considering there is no study on the stability of RG in taro, we supposed that this may be due to actin being a widely used RG.However, it has become evident in some experimental reports that conventional RG, such as actin, is unsuitable owing to its inherent variability [16][17][18].To date, no studies have been published on the validation of RGs for RT-qPCR in taro.
Proper validation of RGs is crucial to ensure the accuracy and reliability of gene expression data, as their expression levels should remain stable under different experimental conditions and across different tissues and growth stages.Therefore, the validation of RGs in taro is a critical step towards understanding its gene expression regulation, identifying key genes associated with important traits, and improving breeding strategies for taro improvement.In this study, we aimed to validate a set of RGs for taro and to provide a reliable toolkit for gene expression studies in this important crop.The findings of this study will contribute to the advancement of taro research and provide valuable insights for the scientific community and stakeholders involved in taro production and utilization (Fig. 1).

Candidate RG screening
To identify suitable RGs for corm growth and different tissues in taro, we employed gene expression value obtained derived from a total of 42 transcriptome datasets from various studies.These datasets encompassed information related to corm color (PRJNA639211) [19], corm growth (PRJNA756190) [20], and four kinds of tissues (Unpublished).The FPKM value (Fragments per kilobase per million mapped reads) of each gene was obtained using StringTie software for gene stable analysis [21,22].To assess the expression stability of each gene, we computed several indices based on the FPKM values using Excel software.These indices included the mean FPKM value (MF), standard deviation (SD), and coefficient of variation (CV) value, which is obtained by dividing the SD by the MF, according to the methods described previously [23].The genes with credible annotation and relative high expression (MF ≥ 10) were selected based on the CV value (≤ 20%).

Plant materials
The taro variety "Lipuyu No. 2" (Colocasia esculenta L. Schott) were planted in field at Guangxi Academy of Agricultural Sciences, Nanning, Guangxi, China.Taro corms at different stages of growth, namely one, two, and three months old, were harvested separately.The leaf, corm, petiole, and root samples were collected at the stage of three months.For drought stress, the samples of three months old were exposed to a seven-days water deficit treatment according to the previous study [24].After harvesting, the samples were thoroughly washed with running water and then with distilled water.Subsequently, the samples were rapidly frozen in liquid nitrogen and stored in a − 80 °C refrigerator until they were processed for RNA extraction.At each sampling time, three samples were collected from three different plants to ensure proper representation for subsequent analysis.

RNA extraction, cDNA synthesis, and qPCR
Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer's instructions.To assess the quality and integrity of the extracted RNAs, agarose electrophoresis and NanoDrop 2000 spectrophotometer (Thermo Scientific, United States) were employed, respectively.To mitigate the risk of DNA contamination, the resulting RNAs underwent treatment with DNase to degrade any residual genomic DNA.Subsequently, one microgram of RNA was utilized for first-strand cDNA synthesis using the HiScript III RT SuperMix for qPCR (Vazyme, China).The ChamQ SYBR qPCR Master Mix (Vazyme, China) was employed for the qPCR reaction.The qRT-PCR assays were conducted in a 25-μL reaction mixture comprising 12.5 μL of Master Mix, 1.0 μL each of forward and reverse primers (0.4 μM), 2 μL of template cDNA (< 100 ng), and 8.5 μL of nuclease-free water.The quantitative amplification was performed using a qTOWER Real-time Thermal Cycler (analytikjena, Germany) with the following cycling conditions: an initial cycle at 95 °C for 30 s, followed by 40 cycles at 95 °C for 5 s, and 60 °C for 30 s.After each amplification, a melting curve analysis was conducted between 55 °C and 95 °C to confirm product specificity.Both no-template and no-RT controls were included for each primer pair.Each RT-qPCR experiment was performed in triplicate.

Analysis of the expression stability of RGs
The expression levels of candidate RGs were assessed based on the cycle of threshold (Ct) value.To determine the expression stability of these candidate RGs, four statistical algorithms were employed: Delta Ct method [25], BestKeeper [26], geNorm [27], and NormFinder [28].For BestKeeper analysis, the raw Ct values were directly utilized, while for the other two analyses, the raw data were transformed into relative quantities using the 2 ^−ΔCt method (where ΔCt = eachCt-minimumCt).To ensure the proper functioning of the four algorithms and to obtain a comprehensive ranking of RGs from Fig. 1 The comprehensive strategy of this study the experimental data, the online tool RefFinder was employed [29].This tool integrated the outputs from the three algorithms, providing a reliable and consolidated assessment of the expression stability of the candidate RGs.

Validation of optimal RGs
To validate the selected RGs, the expression level of a target gene, CeAGPL1 (GeneBank: OK544528), was performed using the same RNA samples employed for RG selection.The CeAGPL1 encode the large subunit of ADP-glucose pyrophosphorylase, which is upregulated in the corm and increases its expression as corm growth progresses [30].cDNA synthesis and qPCR were conducted under the same conditions as described above.
The relative expression level of the CeAGPL1 gene was calculated using the 2^− ΔΔCt method [31].For data normalization, three RG strategies were applied, including: The optimal multiple RGs from each treatment, the most stable RG from each treatment, and the least stable RG.To assess statistical differences, the expression value was subjected to two-way ANOVA using Sidak's test at a significance level of p < 0.05 with a Prism software.For robustness, three replicates of four independent experiments were performed, ensuring reliable and consistent results in the validation process.

Selection of candidate RGs based on the transcriptome data
To select RGs with minimal expression variation for RT-qPCR experiments, we performed a prior analysis of gene stability utilizing RNA-seq data obtained from taro tissues at various stages and treatments.The FPKM values of all transcripts in each taro sample were acquired from transcriptome datasets.The coefficient of variation (CV) value of FPKM for each gene across all samples was calculated and subsequently sorted in ascending order.The gene exhibiting the lowest CV value was chosen as the most promising RG.A total of 14 genes had a CV value lower than 20%.Furthermore, genes should have an average FPKM value higher than 10.Among them, seven annotated genes were selected for analysis of expression stability (Table 1 and Supplementary Table S1).Additionally, three commonly used RGs, namely actin, β-tubulin, and EF-1α, were included.

Specificity and expression analysis of the candidate RGs
To evaluate the specificity and expression profile of the 10 candidate RGs, primers were designed and used to amplify cDNA templates.The melting curves and threshold cycle (Ct) values were determined using RT-qPCR, which represent the specificity and expression level of each gene, respectively.Successful amplification yielded a single band of the expected size, indicating accurate target amplification.Melting curve analysis of primers of all seven candidate genes displayed a single peak, indicating strong primer specificity and the absence of non-specific amplification.The amplification product lengths ranged from 128 to 194 bp.To determine the amplification efficiency of the primers, a serial dilution of cDNA was used as qPCR templates, and the efficiency was calculated based on the slope of the standard curve using Ct values and template concentrations.The calculated efficiencies of the primers ranged from 94.5% to 105.2%.(Table 2, Supplementary Fig. 1).The expression levels of candidate RGs were determined using threshold cycle (Ct) values obtained from RT-qPCR (Supplementary table S2).A smaller Ct value indicated higher gene expression abundance, while a larger Ct value indicated lower abundance.
The average Ct values of each gene ranged from 17.15 (EF-1α) to 27.78 (Actin), which meet the recommended values of higher than 15 and lower than 30 [32].Based on the distribution of raw Ct values, USP25, ACY-1, Actin, CCX4L, and PIA2 exhibited lower variability (Fig. 2A).However, commonly used RGs EF-1α and β-tubulin were considered unfavorable candidates due to their Ct values spanning multiple units, as observed in the interquartile range (Fig. 2B).Nevertheless, relying solely on the distribution of raw Ct values and the interquartile range is insufficient for assessing the expression stability of candidate RGs.Additional methods are necessary for more comprehensive and accurate evaluations.

Expression stability analysis of the candidate RGs
The expression stability of 10 candidate RGs in the leaf, petiole, corm, and root tissues of taro was analyzed under various stages, including different corm growth stages, tissue types, drought stress, and across all samples.The analysis was performed using Delta Ct, BestKeeper, NormFinder, and geNorm methods (Supplementary table S3).

Delta Ct method analysis
Using the Delta Ct method [25], the stability of genes was assessed by calculating the mean standard deviation (SD) of Delta Ct values between one gene and the other 9 genes.A lower mean SD indicates higher gene stability.Table 3

BestKeeper analysis
The BestKeeper method, an Excel-based tool, utilizes pairwise correlations to assess the expression levels of candidate genes [26].We calculated the standard deviation (SD) and coefficient of variance (CV) to evaluate the stability of 10 candidate RGs.The candidate RGs are ranked based on the two indexes with the most stable expression demonstrating the lowest SD value (Fig. 3).The analysis revealed that for corm growth stages, the optimal RG was ACY-1.Actin was identified as the most suitable RG for different tissues, drought stress treatment, and the total sample set.In terms of stability, the least reliable RG was β-tubulin, except for drought stress treatment, in which the most unstable gene was EF-1α.

geNorm analysis
The geNorm program was utilized to calculate the M values of 10 candidate RGs, and the stability of each candidate RG was ranked based on its M value.Lower M values indicate more stable expression, designating the corresponding genes as suitable RGs [27].The samples were categorized into four groups: corm growth stage, different tissues, drought stress, and the overall dataset (all samples).We ranked the candidate RGs based on the average expression stability values (M).The optimal RGs varied depending on experimental conditions.For samples collected at various corm growth stages, ACY-1 and CCX4L were identified as the most reliable RGs (Fig. 4A).USP25 and COX10 exhibited the best stability and were considered optimal RGs among different tissues (Fig. 4B).
In the case of drought stress treatment, PIA2 and Armc8 were deemed the most suitable RGs (Fig. 4C).Lastly, Armc8 and COX10 were identified as the most stable RGs for the entire dataset, encompassing all samples (Fig. 4D).The geNorm algorithm determines the optimal number of RGs by calculating pairwise variation (Vn/n + 1) between consecutive normalization factors, NFn and when Vn/n + 1 is below the threshold of 0.15, the minimum value of n represents the optimal number of RGs required.In this study, Fig. 4E illustrates that the V2/3 values for corm growth stages and different tissues were both lower than the cutoff value of 0.15.This suggests that using two RGs adequately normalized qRT-PCR gene expression data under these two conditions.
However, under drought stress and across all samples, the V3/4 value was lower than 0.15, while the lowest V2/3 value was above 0.15.It suggestes three RGs should be used under these two conditions.

NormFinder analysis
NormFinder is an algorithm that assesses the expression stability of candidate RGs by considering both intra-and inter-treatment expression variations.Lower expression stability values indicate more stable gene expression [28].At corm growth stages and in all samples, PIA2 was identified as the most stable gene, while Armc8 exhibited the highest stability in different tissues and drought stress treatment (Table 4).No wonder, the least stable gene was β-tubulin in all the groups.3 The expression stability of candidate reference genes was assessed using the BestKeeper program.The SD and CV value of each gene were shown in histogram, which represent the expression stability of genes during corm growth stages (A), in various tissues (B), under drought stress (C), and all the samples (D), respectively.SD stands for standard deviation, and CV represents the coefficient of variation.The lower value indicates more stable

Comprehensive ranking of RGs
Upon comparing the outputs of the four methods employed, we noticed some discrepancies in the selection of the most stable RGs.This variance in stability ranking, as observed among Delta Ct, BestKeeper, geNorm, and NormFinder could be attributed to the distinct algorithms utilized by these tools [33,34].To obtain a comprehensive analysis of the results, we calculated the geomean of the ranking values by the RefFinder tool [29].During corm growth, the most stable RGs were ACY-1, expression.An analysis of the pairwise variation (Vn/Vn + 1) between the normalization factors (NFs), denoted as NFn and NFn + 1, were conducted, where n represents the number of genes involved in the NF (E).When the pairwise variation value is lower than 0.15 (dot line), the n number would be regarded as optimal number of genes for normalization Table 4 The RGs were ranked based on their expression stability values determined using the NormFinder analysis PIA2, and CCX4L (Fig. 5A).Among different tissues, the top three most stable RGs were COX10, Armc8, and PIA2 (Fig. 5B).Under drought stress, the most stably expressed RGs were Armc8, COX10, and CCX4L (Fig. 5C).Across all samples, CCX4L, Armc8, and COX10 were the top three most stable RGs (Fig. 5D).Interestingly, the top three RGs for both drought stress and all samples were the same, although the ranking orders were not necessarily identical.By combining the results from geNorm, the optimal number of RGs for normalization was determined.ACY-1 and PIA2 were considered as the optimal multiple RGs for corm growth stages.For different tissues, COX10 and Armc8 were identified as the optimal RGs.However, for drought stress and all the samples, three RGs should be selected for normalization.Specifically, CCX4L, COX10, and Armc8 were utilized as ideal multiple RGs for normalization in these two situations.

Validation of the recommended RGs
To validate the selected candidate RGs, we examined the expression profile of the target gene CeAGPL1, which is involved in increasing starch content.For each treatment, we employed the optimal multiple RGs chosen from that specific condition, as well as the least stable gene, β-tubulin.Additionally, we used the optimal multiple RGs selected from all samples (ACY-1 and PIA2 during corm growth; COX10 and Armc8 within different tissues; and Armc8, COX10, and CCX4L under drought stress) for normalization of the CeAGPL1 expression.The results revealed interesting patterns in CeAGPL1 expression.During corm growth, CeAGPL1 was upregulated continually even based on different RGs (Fig. 6A).However, the expression of the CeAGPL1 gene was significantly upregulated at stage 2 (60 d) using β-tubulin, whereas the changes were not significant using the other two RGs (ACY-1 and PIA2, or only ACY-1).Within different tissues, when normalized with the optimal multiple RGs (COX10 and Armc8) and the most stable gene (COX10), CeAGPL1 expression was significantly increased in corm compared to the leaf (Fig. 6B).However, using the least stable gene (β-tubulin) for normalization between different tissues resulted in downregulation of CeAGPL1, leading to discrepancies compared to normalization with the optimal multiple RGs.Under drought stress, employing the three normalization methods, CeAGPL1 expression exhibited a notable decrease across the three types of tissues (Fig. 6C, D and E).

Discussion
Gene expression analysis plays a crucial role in understanding the growth progress and response mechanisms of plants under various conditions.Several techniques have been developed to investigate gene expression at the transcriptional level, with RT-qPCR being one of the most widely used methods [35].RT-qPCR allows for accurate and reproducible quantification of transcript abundance of target genes.However, it is essential to perform proper normalization using suitable RGs to ensure a reliable assessment of gene expression in RT-qPCR analysis [36].There is mounting evidence indicating that the expression levels of commonly used RGs can exhibit substantial variations in specific contexts.Hence, employing specifically selected RGs that are appropriate for the given conditions is a more reliable approach to obtain accurate gene expression measurements in RT-qPCR [37].With advancements in molecular biology and genomics, taro, a widely grown starchy root crop, has received growing interest in investigating taro at the molecular level [20].Given the absence of validated RGs reported in taro, it becomes crucial to conduct a thorough validation process to identify a set of stable RGs under specific experimental conditions.In this study, RGs suitable for quantifying the expression level of genes involved in corm growth, different tissues, and drought stress were selected and evaluated (Fig. 1).
Screening the most stable expression gene from transcriptome data have been proved an effective strategy [33,34,38].We adopted 42 transcriptome datasets from taro studies that encompassed corm color, corm growth, and tissue.A total of seven potential novel RGs were identified from the transcriptome data of taro based on specific criteria (Table 1 and Supplementary Table S1).Previous studies have demonstrated the reliability of transcriptomic data in identifying suitable RGs in various species, utilizing similar screening criteria [38,39].However, in the case of taro, the stability of actin, which were commonly used as an RG for normalizing RT-qPCR data, has not been evaluated under specific experimental conditions [13][14][15]20].Hence, to assess its expression stability under acid stress, actin was included as a candidate RG in this study.In addition, two other commonly used RGs, EF-1α and β-tubulin [40,41], were employed for validation in these studies.We designed qPCR primers and protocol based on the SYBR dye method with the advantages of ease of use and cost-effectiveness.This choice was made in alignment with our study's principal aim of enhancing accessibility and affordability, ultimately promoting wider participation within taro research community.
The expression levels of these candidate RGs could exhibit significant variations depending on the specific experimental conditions investigated.Among Delta Ct, BestKeeper, geNorm, and NormFinder methods, a slightly divergent pattern in the ranking of expression stability was observed under the same conditions.This finding aligns with previous reports, indicating that variations in results could be attributed to discrepancies in the algorithms and analytical procedures employed by different methods [42][43][44].To counteract potential biases and minimize conflicting outcomes resulting from the use of a single method, we used the RefFinder program to perform a comprehensive ranking of the candidate RGs (Fig. 5).The results revealed that the expression stability of these RGs differed under different conditions, emphasizing the necessity of validating RGs specifically within the context of the experimental treatments employed.Notably, ACY-1, the gene exhibiting the lowest coefficient of variation (CV) value in the transcriptome data, demonstrated remarkable stability across different stages of corm growth.This finding signifies the consistency between validation using qRT-PCR and transcriptome analysis.It has been previously reported that ACY-1 plays a role in regulating root growth and stress response in Nicotiana benthamiana [45].Therefore, it is not surprising that ACY-1 displayed stable expression during corm growth and exhibited relatively lower stability under drought stress conditions in our study.PIA2 and CCX4L were ranked as the second-and third-most stable genes, respectively, during corm growth.Interestingly, COX10 and Armc8 were identified as the top three stable genes in the other three groups, including different tissues, drought stress, and all samples.Moreover, the three most stable genes under drought stress conditions (Armc8, COX10, and CCX4L) were consistent with the top three stable genes observed across all the samples.However, there was no single universal RG that consistently ranked among the top three across all experimental groups.The expression stability of RGs varied depending on the specific experimental conditions and treatments being studied.This highlights the importance of selecting appropriate RGs and validating their stability for each experimental scenario to ensure accurate and reliable normalization of gene expression data.
In addition to evaluating the expression stability of candidate RGs and identifying the most stable ones, another crucial aspect of RG validation is determining the optimal number of RGs for proper normalization [46].Studies have highlighted the necessity of using multiple RGs for accurate assessment of gene expression under specific experimental conditions [47][48][49].The rationale behind using multiple RGs lies in the assumption that any experimental errors present would affect all genes expressed in the same sample, assuming they are processed together.Consequently, the individual replicates' experimental errors are averaged across the RGs, and employing the geometric mean of multiple RGs yields a more robust estimate of the overall experimental error compared to using individual RG.The geNorm algorithm is frequently employed to ascertain the optimal number of RGs through analysis of pairwise variation (Vn/n + 1) [27].For corm growth and different tissue groups, the V2/3 value was below the cut-off of 1.5, indicating that two RGs would suffice to normalize gene expression.Therefore, the recommended RGs for corm growth and different tissues were ACY-1/PIA2 and COX10/Armc8, respectively (Fig. 5 A and B).However, under drought stress, the V4/5 value was less than 0.15, suggesting the need for at least four RGs to accurately evaluate gene expression.In all the sample groups, the first value that lower than 0.15 was V3/4, indicating that using of three RGs is more accurate (Fig. 4).Given the increased cost and effort associated with using a higher number of RGs, it is reasonable to consider utilizing three RGs (Armc8, COX10, and CCX4L) that demonstrated stability for both drought stress and all samples (Fig. 5C and D).
In our study, three commonly used RGs were employed for stability analysis.However, EF1_α and β-tubulin were ranked as the least stably expressed genes (Fig. 5).Several recent studies have highlighted that the expression levels of certain traditional RGs can exhibit significant variations under specific experimental conditions, challenging the assumption that they are inherently stable, as previously believed [50,51].These findings emphasize the importance of critically re-evaluating the suitability of RGs and highlight the need for careful selection and validation of appropriate RGs tailored to the specific experimental context for accurate and reliable normalization of gene expression data.
To validate the selected RGs, three RG strategies were applied to normalize the expression level of the target gene CeAGPL1, which encodes an AGPase large subunit involved in increasing starch content [30,52].When normalized using the optimal RGs, except under drought stress, the expression of was notably upregulated in corm growth stage 3, but not in stage 2 (Fig. 6A).In addition, with optimal RGs, the transcriptional level of CeAGPL1 was significantly higher in corm than in leaves (Fig. 6B).Similar expression trends of CeAGPL1 were observed in our previous study, where the expression of CeAGPL1 was accessed by transcriptome [30].However, when normalized with the least stable RG of β-tubulin, overestimation of CeAGPL1 expression was observed in stage 2, and the expression of CeAGPL1 in corm was underestimated, although not significantly.The choice of RGs significantly influenced the expression pattern of CeAGPL1.This highlights the critical importance of using reliable RGs as a prerequisite for accurate RT-qPCR data analysis of taro under various conditions.Furthermore, it is important to acknowledge the limitations of the methods used for RNA quantification in this study.Although we primarily employed a spectrophotometer for RNA quantification, it is essential to recognize that this method may have shortcomings in accurately quantifying RNA.Alternative methods, such as a fluorimeter with RNA-specific dyes, would offer greater accuracy in RNA quantification, and future studies may benefit from incorporating these methods to enhance the accuracy of RNA quantification.

Conclusion
This study highlights the potential of transcriptomic datasets from taro to identify novel candidate RGs for RT-qPCR assays and emphasizes the importance of experimental validation of gene expression stability at diverse stages.Among seven candidate genes with relatively high expression levels and low variance, transcriptomic datasets generated during corm growth and different taro tissues were used for screening.Three commonly used RGs were also included in validation.Through the application of four statistical algorithms (Delta Ct, BestKeeper, geNorm, and NormFinder) and a comprehensive program (RefFinder) for analysis, along with validation using the target gene ACY-1 and PIA2 were identified as ideal RGs during corm growth, while COX10 and Armc8 were suitable for different tissues, Armc8, COX10, and CCX4L were the optimal RGs under drought stress.Significantly, it was observed that β-tubulin was not stably expressed in taro.This study represents the first systematic analysis of RGs for RT-qPCR in taro, especially during corm growth.These findings are valuable for future investigations into functional genes related to corm growth in taro and will aid in selecting suitable RGs for starch biosynthesis research in tuberous crop plants during starch accumulation.

u b u l i n 15 Fig. 2
Fig. 2 RT-qPCR Ct values and interquartile ranges are shown.A Ct values of each RG in all samples are presented.The median is represented by a line across the box, while the box indicates the 25th and 75th percentiles.Whiskers represent the maximum and minimum values.B The interquartile ranges demonstrate the variability of Ct values between the 25th and 75th percentiles Fig.3The expression stability of candidate reference genes was assessed using the BestKeeper program.The SD and CV value of each gene were shown in histogram, which represent the expression stability of genes during corm growth stages (A), in various tissues (B), under drought stress (C), and all the samples (D), respectively.SD stands for standard deviation, and CV represents the coefficient of variation.The lower value indicates more stable

Fig. 4
Fig.4 The stability and optimal number of RGs evaluated by geNorm software.The stability ranking of each group, as shown in A (Corm growth), B (Different tissues), C (Drought stress), and D (All samples), were presented in line chart.A lower value of the M value indicates more stable expression.An analysis of the pairwise variation (Vn/Vn + 1) between the normalization factors (NFs), denoted as NFn and NFn + 1, were conducted, where n represents the number of genes involved in the NF (E).When the pairwise variation value is lower than 0.15 (dot line), the n number would be regarded as optimal number of genes for normalization

Fig. 5
Fig.5 The comprehensive ranking of RGs in four groups using RefFinder program.The four groups are A (Corm growth), B (Different tissues), C (Drought stress), and D (All samples).The gene with lower Geomean value is more stable

S t a g e 1 SFig. 6
Fig.6 Relative expression of CeAGPL1 normalized with different sets of RGs.The relative expression of CeAGPL1 was validated in three corm growth stages (A), across four kinds of tissues (B), and under drought stress in corm (C), leaf (D), and root (E).Three normalization strategies were used, including multiple top-ranked RGs (Blue color), single top-ranked RG (Green), and the least stable RG (Red).The expression change with significant was indicated in line and asterisk

Table 1
The summary information of candidate RGs in taro based on the transcriptome data

Table 2
Primer sequences of 10 genes used for RG validation

Table 3
The RGs were ranked based on their suitability for normalization, and their expression stability values were determined using the Delta Ct method